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This paper presents novel extensions and applications of the Uppaal-SMC model checker. The 
extensions allow for statistical model checking of stochastic hybrid systems. We show how our 
race-based stochastic semantics extends to networks of hybrid systems, and indicate the integration 
technique applied for implementing this semantics in the Uppaal-SMC simulation engine. We re- 
port on two applications of the resulting tool-set coming from systems biology and energy aware 
buildings. 

1 Introduction 

Statistical Model Checking (SMC) |l29l|27l[32l|34l|25l is an approacli that has recently been proposed as 
new validation technique for large-scale, complex systems. The core idea of SMC is to conduct some 
simulations of the system, monitor them, and then use statistical methods (including sequential hypoth- 
esis testing or Monte Carlo simulation) in order to decide with some degree of confidence whether the 
system satisfies the property or not. By nature, SMC is a compromise between testing and classical for- 
mal method techniques. Simulation-based methods are known to be far less memory and time intensive 
than exhaustive ones, and are some times the only option. SMC has been implemented in a series of tools 
that have defeated well-known tools such as PRISM on several case studies. Unlike more "academic" 
exhaustive (and intractable) techniques, SMC is spreading to various research areas such as systems 
biology |[2Tll24l and software engineering |[36l[30l . in particular for industrial applications HID HI. 

There are several reasons for this success. Firstly, SMC is simple to understand, use and (in princi- 
ple) to implement. Secondly, no additional modelling or specification effort is needed, provided that the 
mode-ling formalism used can be given a natural stochastic semantics serving as the basis for interpre- 
tation of the specification formalism and as a basis for generating simulation-runs. Thirdly, SMC allows 
to analyse properties 1 1 1 , 4 | that cannot be expressed in classical temporal logic, including properties for 
which classical model checking is undecidable. 

In a series of recent works |[T6l[T5l . we have investigated the problem of Statistical Model Checking 
for networks of Priced Timed Automata (PTAs). PTAs are timed automata, whose clocks can evolve 
with different rates, while0 being used with no restrictions in guards and invariants. In |T5l, we have 
proposed a natural stochastic semantics for such automata, which allows to perform statistical model 

* Work partially supported by the VKR Centre of Excellence MT-LAB, the Sino-Danish Basic Research Center IDEA4CPS, 
and by the CREATIVE project ESTASE. 

'in contrast to the usual restriction of priced timed automata fTllTl 



Ezio Bartocci and Luca Bortolussi (Eds.): HSB 2012 
EPTCS 92, 2012, pp. 1224T36] doi: 10.4204/EPTCS .9Z91 



© David, Du, Larsen, Legay, Mikucionis, Poulsen, Sedwards 
This work is licensed under the 
[ Creative Com mons Attrib iition License. 



David, Du, Larsen, Legay, Mikucionis, Poulsen, Sedwards 



123 



checking. Our work has later been implemented in Uppaal-SMC, that is a stochastic and statistical 
model checking extension of Uppaal. Uppaal-SMC relies on a series of extensions of the statistical 
model checking approach generalised to handle real-time systems and estimate undecidable problems. 
Uppaal-SMC comes together with a friendly user interface that allows a user to specify complex prob- 
lems in an efficient manner as well as to get feedback in the form of probability distributions and compare 
probabilities to analyse performance aspects of systems. The Uppaal-SMC model checking has been 
applied to a wide range of examples from networking and Nash equilibrium! 9] through systems biol- 
ogy l[T4l . real-time scheduling |fT3l . and energy aware systems |[T2]| . 

For PTAs clocks evolve with fixed rates depending only on the discrete state of the model, i.e. loca- 
tions and discrete variables. In this paper, we present and implement an extension of the modelling for- 
malism of Uppaal-SMC, where clock rates may depend not only on values of discrete variables but also 
on the value of other clocks, effectively amounting to ordinary differential equations (ODEs). As a first 
contribution of this paper, we present an extension of our race-based stochastic semantics to networks of 
stochastic hybrid automata (SHA). One major difficulty is the implementation of this extensions in our 
engine for generating random runs, in a manner which is correct with respect to the stochastic semantics. 
In fact Uppaal-SMC does not solve the equations exactly but currently supports the Euler integration 
method. A fixed time step (defined by the user) is used by an internal integrator component added to the 
system. This integrator races with the other processes during its time step and all rates ai^e considered 
constant and defined by the equations in the model. As a second contribution of the paper we describe 
this integration in Uppaal-SMC. It is worth mentioning that while other SMC -based approaches exist 
for SHAs |[36l[35l l2]|. none of them do consider ODE-based modelisation for clock rates. 

To demonstrate the applicability of this new extension, we apply advanced SMC techniques to two 
challenging applications coming from systems biology and energy aware buildings. For the case of 
systems biology, we show how the combination of ODEs and SMC allows us to reason on biological os- 
cillations - a problem that is beyond the scope of most existing formal verification techniques. We model 
a genetic circadian oscillator, which is used to distil the essence of several real circadian oscillators. For 
the case of energy aware buildings, we refer to a recently developed framework including components 
for layout of buildings, availability of heaters, climate and user behaviours allowing to evaluate different 
strategies for distributing heaters among rooms in terms of the resulting comfort and energy consumption. 
To indicate central parts of this framework and the clear advantages of modelling the evaluation of room 
temperatures with ODEs, we illustrate in this paper the framework with a small instance comprising two 
rooms with a single shared heater. 

Structure of the Paper. The remainder of the paper is structured as follows. In the next Section |2] 
we preview the expressive power of the hybrid extensions of UPPAAL-SMC using an extension of the 
well-known bouncing ball. Sections [3] and |4] details the semantics of the extended formalism of networks 
of stochastic hybrid automata. Applications to Energy Aware Buildings and a Biological Oscillator are 
given in Section [5] and Section [6l respectively. Finally, Section |7] concludes the paper and suggests 
directions for future research. 



2 Throwing, Bouncing and Hitting Ball 

To give an early illustration of the expressive power of the extended modelling formalism of Uppaal- 
SMC, we consider a variant of the well-known bouncing ball. In our versioiH the ball is initially thrown 
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x = 10, 

vx = -10+random(0.5) 

bounce 
X == && vx < 
vx = -(0.85+random(0.15))*vx 

(a) Model of x coordinate. 
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vy' == -9.81 &&. 
y=10, y'==1*vy 
vy = 5.8+random(0.5) 



bounce! 




y == && vy < 
vy = -(0.8+random(0.12))*vy 



hit? 

y >= 6 && vy >= && 
x>= 10 

vy=-(0.85+random(0.1))*vy-4 
hit? 

y>=6 && vy<0 && vy>=-4 && 
x>=10 



vy = -4.0 
(b) Model of y coordinate. 



hit! 



(c) Player trying to hit with exponentially 
distributed delay. 




Qy 



(d) Trajectory as {x,y) plot. 
Figure 1: Models and a trajectory of a thrown/bouncing ball hit by a player. 



against a wall, bounces against it, and then continues its trajectory by falling and bouncing against the 
floor. In addition, a (inexperienced) player tries to hit the ball randomly according to an exponential 
distribution. The model is depicted in Fig. [Ila)-(c). The player is modelled as a simple automaton that 
broadcasts hit ! with an exponential distribution of rate 5/2. The x coordinate (Fig. [Tall is initialised to 
10 with an uncertain derivative vx uniformly between [—10, —9.5] (the ball is thrown against the wall), 
after which the ball moves toward the wall (placed at 0). Here, the automaton outputs bounce ! on an 
urgent channel, which forces the transition to take place deterministically at x= 0. After a bounce with 
a random dampening factor of the velocity vx uniformly between [0.85, 1], the ball continues to move in 
the opposite direction. The y coordinate (Fig.fTbl) is initialised to 10 with an uncertain derivative vy. The 
model shows the effect of gravitation with vy'= —9.81. The ball bounces with a random dampening 
factor on the floor (at 0) and when the ball is away from the wall (x> 10) then it can be hit by the player 
provided it is high enough (y> 6). Depending on the current direction of the ball, the ball may bounce or 
it is pushed. One possible trajectory of the ball is shown in Figure [Td] The plot is obtained by checking 
the query "simulate 1 [x<=40] {y}". The vertical line shows the ball moving to its initial position 
and should be ignored. The ball bounces as expected against the wall, the floor, and the hitting of the 
player. Uppaal-SMC is able to simulate this hybrid system that has a second order ODE, a stochastic 
controller (the player), and a stochastic environment (random dampening factor). 

In addition, we may perform statistical model-checking in order to estimate the probability that the 
ball is still bouncing above a height of 4 after 12 time units with the query: 

Pr[<=20](<> tinie>=12 and y>=4) 

which returns the confidence interval [0.44,0.55] with 95% confidence after having generated 738 runs. 
We can also test for the hypothesis 



^Our version is inspired by the recently announced solution by a 16 year old German school boy (Shouryya Ray), to a 350 
year open problem by Newton concerned with predicting the trajectory of a ball thrown at a wall. 
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Pr[<=20](<> tinie>=12 and y>=4) >= 0.45, 

which gives a more precise lower bound. The hypothesis holds with a region of indifference ±0.01 and 
a level of significance of 5% after generating 970 runs. 

3 Networks of Hybrid Automata 

Uppaal-SMC now supports the analysis of stochastic hybrid automata (SHA) that are timed automata 
whose clock rates can be changed to be constants or expressions depending on other clocks, effectively 
defining ODEs. This generalizes the model used in our previous work |[T6l[T5l where only linear priced 
automata were handled. Our new release Uppaal-SMC 4. l.lcH supports fully hybrid automata with 
ODEs and a few built-in complex functions (such as sin, cos, log, exp and sqrt)! 

Hybrid Automata Intuitively, a hybrid automaton is a finite-state automaton extended with con- 
tinuous variables that evolve according to dynamics characterizing each discrete state (called a location). 
Let X be a finite set of continuous variables. A variable valuation over Z is a mapping v : X ^ M, where 
M is the set of reals. We write for the set of valuations over X. Valuations over X evolve over time 
according to delay functions F : ]R>o x — )• M^, where for a delay d and valuation v, F{d, v) provides 
the new valuation after a delay of d. As is the case for delays in timed automata, delay functions are as- 
sumed to be time additive in the sense that F{d\,F{d2, v)) = F(<ii v). To allow for communication 
between different hybrid automata we assume a set of actions £, which is partitioned into disjoint sets of 
input and output actions, i.e. £ = Z,- l±l Lo- 

Definition 1. A Hybrid Automaton (HA) .^^ is a tuple = {L,Iq,X ,E ,F,I), where: (i) L is a finite set 
of locations, ( ii) Iq £L is an initial location, ( Hi) X is a finite set of continuous variables, (iv)'L = 1+) 
is a finite set of actions partitioned into inputs and outputs (T,o), (v) E is a finite set of edges of 
the form {£,g,a,(p,£'), where £ and £' are locations, g is a predicate on R^, action label a € T and (p 
is a binary relation on M.'^, (vi)for each location £ £ L F{£) is a delay function, and (vii) I assigns an 
invariant predicate I{£) to any location £. 

The semantics of a HA M' is a timed labeled transition system, whose states are pairs {£, v) € L x 

with V \= I{£), and whose transitions are either delay transitions {£,v) i^,v') with d G R>q and 
v' = F{d,v), or discrete transitions {£,v) {£',v') if there is an edge {£,g,a,(p,£') such that V \= g 
and 9(v, v'). We write {£, v) v') if there is a finite sequence of delay and discrete transitions from 

(Av)to(f,v'). 

In the above definition, we have deliberately left open the concrete syntax for the delay function 7^ 
as well as the invariant I. In Uppaal-SMC the delay update for a simple clock x - used in (priced) timed 
automata - is given by an implicit rate x' = 1 or an explicit rate x' = e appearing in the invariant of £, 
where e is an expression only depending on the discrete part of the current state. More generally, the 
effect of the delay function F may be specified by a set of ODEs needing to be solved. It is important to 
note that in specifying the delay function F and the invariant 7, the full syntax of Uppaal expressions - 
including user-defined functions - is at the disposal. 

Example 1. Reconsider the extended bouncing ball example from Section |2] Here the automaton for 
the X coordinate may be initialized to the state (x = 10, vx = —9.8), after which the following transition 
sequence may occur: 

(x = 10, vx = -9.8) (x = 0,vx = -9.8) (x = 0,vx = 9.31) 

■'www . uppaal . org. 
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where in the bounce! -transition the dampening factor has non-deterministically been chosen from the 
interval [0.85, 1.00] as 0.95. The automaton for the y coordinate may be initialized to the state (y = 
10, vy = 5.9) after which the delay function will effectively be given by delay transitions: 

(y= 10,vy = 5.9) A (y = -9.81/2j2 + 5.9J+10,vy = -9.8W + 5.9) 

The ball bouncel'es on the floor, i.e. y = or when (approximately) d = 2.15, and will afterwards 
- having non-deterministically chosen the dampening factor 0.9 - continue from the state (y = 0,vy = 
13.67). 

Networks of Hybrid Automata Following the compositional specification theory for timed systems in 
ifTTl . we shall assume that NHAs are input-enabled in the sense, that for all states {I, v) and input actions 
I G for all HAs j, there is an edge l, <p,£^') such that V \= g and <p(v, v') for some valuation v'. 
Also, we shall assume that time always diverges, and that different automata synchronize on matching 
inputs and outputs as a standard broadcast synchronization [20 1 ■ 

Whenever jz/^ = {U ,Xj ,Lj ,Ej ,Fj jj) (j = l...n) are NHA, they are composable into a closed 
network iff their variable-sets are disjoint {X^ n X*^ = when j ^ k), they have the same action set 
(E = E-' = for all j,k), and their output action-sets provide a partition of £ (Li n = for j ¥^ ^» ^rid 
E = UjZo)- For a € r we denote by c(a) the unique j with a G If V G with X = UjX^, we denote 
by V ixj the projection of V to X->. 

Definition 2. Let £/j = ,Xj ,F-> ,!■>) (with j = I . . .n) be composable NHAs. Their composition 

I • • • I is the HA£/ = {L,X,'L,E,F,I) where (i) L = XjV, (ii) X = UjX^, (Hi) F{l){d, v){x) = 
Fj(^ij){d,v ixj)i^) whenxeXj, (iv) = njl{i^), and(v) {i,njgj,a,Ujq)j,i') eE whenever {ij,gj,a,q)j, 
EJ for j = \ ...n. 

4 Stochastic Semantics for Networks of Hybrid Automata 

Reconsidering again our extended version of the bouncing ball from Section [2l it is clear that there is 
a constant race between the ball bounce ling on the floor and the player hit ling the ball. Whereas 
the time of bouncing is deterministic - given by the ODE obtained from the (stochastic) effect of the 
previous bounce I or hit I - the time of hitting is stochastic according to an exponential distribution 
with rate 5/2. However based on this, a measure on sets of runs of the systems is induced, according 
to which quantitative properties such as "the probability that the ball with have a height greater than 4 
after 12 time-units" become well-defined. 

Our early works |[T5l - though aimed at stochastic semantics of priced timed automata - is suffi- 
ciently general that it also provides the basis for a natural stochastic semantics for networks of HAs, 
where components associate probability distributions to both the time-delays spent in a given state as 
well as to the transition between states. 

Let si^ = [L-i ,X-i ,L,E-i ,F^ ,H) (7 = 1 . . . «) be a collection of composable HAs. Under the assump- 
tion of input-enabledness, disjointedness of clock sets and output actions, states of the composite NHA 
j2/ = I ... I may be seen as tuples s = (si,... ,Sn) where sj is a state of £/ i.e. of the form {£, v) 
where i G and v G M^^ . Our probabilistic semantics is based on the principle of independence between 
components. Repeatedly each component decides on its own - based on a given delay density function 
and output probability function - how much to delay before outputting and what output to broadcast at 
that moment. Obviously, in such a race between components the outcome will be determined by the 
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component that has chosen to output after the minimum delay: the output is broadcast and all other 
components may consequently change state. 

Stochastic Semantics of HA Components The stochastic semantics of HAs refine the non-deterministic 
choices that may exist with respect to delay, output and next state. We consider the component £/ ^ and 
let St^ denote the corresponding set of states. For each state s = {£, v) of £/■', we shall assume that there 
exist probability distributions for delays, output as well as next-state: 

• the delay density function, /i, over delays in M>o, provides stochastic information for when the 
component will perform an output, thus / lJLs{t)dt = 1; 

• the output probability function 7^ assigns probabilities for resolving what output o € S/j to generate, 
i-e- Lo7.(o) = 1; 

• the next-state density function r]f provide stochastic information on the next state s' = {i' , V) G M'''^ 
given an action a, i.e. f^i ri"{s') = 1. 

For outputs happening deterministically at an exact time point d (or deterministic next states s'), /i^ (rj") 
becomes a Dirac delta function 5j (4' D- 

In Uppaal-SMC uniform distributions are applied for states where delay is bounded, and exponen- 
tial distributions (with location-specified rates) are applied for the cases, where a component can remain 
indefinitely in a location. Also, Uppaal-SMC provides syntax for assigning discrete probabilities to dif- 
ferent outputs as well as specifying stochastic distributions on next-states (using the function random [b] 
denoting a uniform distribution on [0,^]). 

Stochastic Semantics of Networks of HA For the stochastic semantics of closed networks of HA con- 
sider £/ = I ... I £/n) with a state space St = Sti x • • • x St„. For s = (i'l , . . . G St and aia2 ...ak& 
£* we denote by n{s,a\a2 . ..au) the set of all maximal runs from s with a prefix tia\t2a2 ■ ■ ■ tuau for some 
ti,... ,tn G M>o, that is runs where the i'th action a,- has been outputted by the component Ac(a,). Pro- 
viding the basic elements of a Sigma-algebra, we now inductively define the following measure for such 
sets of runs: 

P.c/(7r(s,ai . ..«„)) = 

/ M...(o-(n/ /(n<;'(-o-p^(^(s':«2...fl«))«'s')«'f 

Jt>0 j^^JT>t Js' j J 

This definition requires a few words of explanation: at the outermost level we integrate over all 
possible initial delays t. For a given delay t, the outputting component c = c{ai) will choose to make 
the broadcast at time t with the stated density. Independently, the other components will choose to a 
delay amount, which - in order for c to be the winner - must be larger than t; hence the product of the 
probabilities that they each make such a choice. Having decided for making the broadcast at time t, the 
probability of actually outputting ai is included. Finally, integrating over all global states s' that may 
result from all components having delayed t time-units and changed state stochastically with respect to 
the broadcasted action ai , the probability of runs according to the remaining actions 02 • ■ • is taken into 
account. 



which should formally be treated as the limit of a sequence of delay density functions with decreasing, non-zero support 
around d. 
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Generation of Random Runs The central component of the SMC engine of Uppaal-smc is the 
efficient generation of random runs according to the stochastic semantics proposed in the preceding 
Section. Uppaal-SMC has to integrate the ODEs given in the model while still respecting its stochastic 
semantics. The tool does not solve the equations exactly and currently supports the Euler integration 
method. A fixed time step 5t (defined by the user) is used by an internal integrator component added to 
the system. This integrator races with the other processes during its time step and all rates are considered 
constant and defined by the equations in the model. At the end of every step or if another process wins 
the race, all the rates are re-evaluated and the resulting values of the clocks are computed according to 
Xnew = ^oid + St *x'gi^j. We plan to implement more robust methods, such as Runge-Kutta's, that give 
different evaluations of x„e».. In this case, x'^^ is then derived from the method and not directly from the 
equation in the model and kept constant during 5t. The stochastic semantics and the races between the 
components are still respected. 

The integrator component always picks a delay equal to the given time step St. The other components 
pick delays according to their distributions. The winning component delays to the point in time it wanted 
and tries to take a transition at that point. If it appears that no guard is satisfied, no transition is taken 
and all components race again. In practice, every component will use the current rates of the clocks to 
evaluate a lower bound on the delay from which anything can happen to skip delaying for nothing. 

Statistical Model Checking We use SMC (HI [32l |34l |29l to estimate and test on the probability that 
a random run of a network of SHAs will satisfy a given property. Given a model and a trace property 
(p (e.g. expressed in LTL |[3T| or MTL |[26l ). SMC refers to a series of simulation-based techniques 
that can be used to answer two questions: (1) Qualitative: is the probability that a random run of 
will satisfy (p greater or equal to a certain threshold 6 (or greater or equal to the probability to satisfy 
another property (p')7 and (2) Quantitative: what is the probability that a random run of will satisfy 
<p? In both cases, the answer will be correct up to a user-specified level of confidence, providing a upper 
bound on the probability that the conclusion made by the algorithm will be wrong. For the quantitative 
approach, which we will use intensively in this paper, the method computes a confidence interval that 
is an interval of probabilities that contains the true probability to satisfy the property. Here the confi- 
dence level provides the probability that the computed confidence interval indeed contains the unknown 
probability. 

Our Uppaal-SMC tool-set implements a wide range of SMC algorithms for networks of SHAs not 
only for reachability and safety properties, but also for general weighted MTL properties ifTOl [8l. In 
addition, the tool offers several features to visualize and reason on the results. 

5 Energy Aware Buildings 

Uppaal-SMC has recently (12] been applied to an evaluation framework for energy aware buildings, 
where the control of heating is optimized with respect to environmental and user profiles. To indicate 
central parts of this framework and the clear benefit of modeling with ODEs, we illustrate the framework 
on a simplified instance and recall results from our case study lfT2l based on a benchmark for hybrid 
systems verification lilSJ . 

A Simple 2-Room Example To illustrate the various aspects of the (extended) modeling formalism 
supported by Uppaal-SMC, we consider the case of two independent rooms that can be heated by a 
single heater shared by the two rooms, i.e., at most one room can be heated at a time. Figj^a) shows 
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0N_1 

(a) stochastic heater. 
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T[0]=TO[0; 
OFF 



on[0]? 

K =9+rand om(3) 



. offfOI? 



Init 



ON 



T[l]=TO[i: 
OFF 



on[l]? 

K =9+rand om(3) 



. off [11? 



ON 



T[0]" = = -T[0]/10 K" = = 0&& T[l]" = = -T[l]/10 K" ==0&& 
+sumy:int[0,l]) T[01"= = K - T[0]/10 +sumg:int[0,l]) T[1]"= = K - T[l]/10 

A[0]|j]*(TG]+-T[0])/s +sum(j:int[0,l]) A[l]D]*(T[j]+-T[l])/s +sum(j:int[0,l]) 

A[0][j]*(T[j]+-T[0])/s A[l][j]*(TU] + -T[l])/s 



(b) room 0. 



(c) room 1. 



Figure 2: A simple two room example with an autonomous heater. 



the automaton for the heater. It turns itself on with a uniform distribution over time in-between [0,4] 
time units. With probability 1/4 room is chosen and with probability 3/4 room 1. The heater stays 
on for some time given by an exponential distribution (rate 2 for room 0, rate 1 for the room 1). In 
summary, one may say that the controller is more eager to initiate the heating of room 1 than room 
0, as well as less eager to stop heating room 1. The rooms are similar and are modeled by the same 
template instantiated twice as shown in Fig.[3b-c). The room is initialized to its initial temperature and 
then depending on whether the heater is turned on or not, the evolution of the temperature is given by 
T/ = -Ti/lO + j:j=o,iAij{Tj - Ti) or T! = K - 7]/10 + Zj=Q^iAij{Tj - 7]) where ij = 0, 1 are room 
identifiers. The sum expression corresponds to an energy flow between rooms and matrix A encodes 
the energy transfer coefficient between adjacent rooms. Furthermore, when the heater is turned on, its 
heating is not exact and is picked with a uniform distribution of K € [9,12], realized by the update 
K=9+random(3). 

This example illustrates the support for NSHA in Uppaal-SMC with extended arithmetics on clocks 
and generalized clock rates. 



Extended Input Language Uppaal-SMC takes as input NSHA as described above. Additionally, 
there is support for other features of the Uppaal model checker's input language such as integer vari- 
ables, data structures and user-defined functions, which greatly ease modeling. Uppaal-SMC allows the 
user to specify arbitrary rates for the clocks, which includes a mix of integer and clock expressions on any 
location. In addition, the automata support branching edges where weights can be added to give a distri- 
bution on discrete transitions. It is important to note that rates and weights may be general expressions 
that depend on the states and not just simple constants. 



Checking Queries The fundamental principle in Uppaal-SMC is to generate runs and evaluate some 
expression on the states along the obtained run. Runs are always bounded, either by time, by a number 
of steps, or more generally by cost (when using a clock explicitly). The engine has a built-in heuristic 
detection of Zeno behaviours to abort the generation of such runs. Examples of the syntax for the 
different types of bounds are [<=100] for 100 time units since the beginning of the run, [#<=50] for 50 
discrete transitions taken from the initial state, and [x<=200] until the clock x reaches 

Uppaal-SMC supports simulations with monitoring custom expressions, probability evaluation, hy- 
pothesis testing, and probability comparison. We can simulate and plot the temperatures with the query 

simulate 1 [<=600] {T [0] ,T [1] } 



^It is up to the modeler to ensure that the clock eventually reaches the bound. 
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The query request the checker to provide one simulate run over 600 time units and plot the temperatures 
of Room(O) and Room(l). The heater in this example is purely stochastic and is not intended to enforce 
any particular property. Yet, the simulation obtained from this query in Fig. [3] shows that the heater is 
able to maintain the temperatures within (mostly) distinct intervals. 

We can evaluate on a shorter time scale the probability for the temperature of Room (0) to stay below 
30 and the temperature of Room(l) to stay above 5 with the queries 

Pr[<=100]([] Room(O) .Init || T[0] <= 20) 
Pr[<=100]([] Room(l) .Init || T[l] >= 7) 

The results are respectively in [0.45,0.55] and [0.65,0.75]. The precision and confidence of these so- 
called confidence intervals are user-defined and influence the number of runs needed to compute the 
probability. In this example, for having the precision to be ±0.05 with a confidence of 95%, we need 
738 runs. In fact if we are only interested in knowing if the second probability is above a threshold it 
may be more efficient to test the hypothesis 

Pr[<=100]([] Room(l) .Init || T[l] >= 7) >= 0.69 

which is accepted in our case with 902 runs for a level of significance of 95%. To obtain an answer at 
comparable level of precision with probability evaluation, we would need to use a precision of ±0.005, 
which would require 73778 runs instead. 

The tool can also compare probabilities without needing to compute them individually. We can test 
the hypothesis that the heater is better at keeping the temperature of Room(l) above 8 than keeping the 
temperature of Room (0) below 20: 

Pr[<=100]([] Room(l).Init || T[l] >= 7) >= 
Pr[<=100]([] Room(0).lnit || T[0] <= 20) 

which is accepted in this case with 95% level of significance with just 258 runs. 



Results In |[T2l we have estimated the comfort time (duration while being in comfortable temperature 
range) and energy consumption for various weather conditions, user profiles and central controller strate- 
gies. Fig.|4]shows six energy consumption estimates in different configurations (a building with 5 rooms 
and 3 heaters). The energy comparison shows that the dynamic user profile can save more than 33% of 
energy regardless of the chosen central controller strategy. 
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Figure 4: Energy consumption estimates for static and dynamic user profiles. 



6 Biological Oscillator 



One of the key oscillatory behaviours in biology is the circadian rhythm that allows an organism to 
take advantage of periods of day and night to optimise when to maximise activity and recovery. We 
show how the genetic circadian oscillator of |l3l|33l can be modelled and analysed using Uppaal-SMC. 
This synthetic model distils the essence of several real circadian oscillators and demonstrates how a 
reliable system can be constructed in the face of inherent stochasticity. Figure |5a] shows a system of 
differential equations from |[33]| . The equations are typeset in Uppaal as invariant expression on a 
location shown in Fig. [5bJ where each quantity DA,Dii,D'^,MA,MR,A,R,C are modelled as continuous 
clock variables DA, DR, D_A, MA etc with rates defined by a corresponding differential equation. The 
preceding expressions about alphaA, alpha_A, alphaR, alpha_R, beta_A and so on are modelling the 
constants o^, a^, a^, a^, jSa and so on. The assignments on the first transition initialise all the variables 
with initial conditions. Uppaal-SMC is then used to simulate the model and provide plots of how the 
variable values evolve over time, which are displayed in Figure ITal 



dDA/dt 


= QaD^a - JaDaA 


dDR/dt 


= QrD'r - YrDrA 


dD\ldt 




dD'^/dt 


= rRDRA-0RD'R 




= a'^D'^ + UADA-SMf.MA 


dMR/dt 


= a'RD'R + UrDr - 5msMr 


dA/dt 


= PaMa + BaD'a + QrD'r 




- A(yADA + yRDR + ycR + SA) 


dR/dt 


= PrMr-YcAR + SaC-SrR 


dC/dt 


= YcAR-SaC 



(a) Ordinary differential equations. 



alphaA=50, alpha_A=500, alphaR=0.01, alpha_R=50, 
betaA=50, betaR=5, deltaMA=10, deltaMR=0.5, deltaA=1, deltaR=0.2, 
gammaA=1, gammaR=1, gammaC=2, thetaA=50, thetaR=100, 
< DA=1, DR=1, D_A=0, D_R=0, MA=0, MR=0, A=0, R=0, C=0 

Cj alphaA'==0 && alpha_A'==0 && alphaR'==0 && alpha_R'==0 && 
betaA'==0 && betaR'==0 && deltaA'==0 && deltaR'==0 && 
deltaMA'==0 && deltaMR'==0 && gammaA'==0 && 
gammaR'==0 && gammaC'==0 && thetaA'==0 && thetaR'==0 && 
DA'== thetaA*D_A-gammaA*DA*A && 
DR'== thetaR*D_R-gammaR*DR*A && 
D_A'== gammaA*DA*A-thetaA*D_A && 
D_R'== gammaR*DR*A-thetaR*D_R && 
MA'== alpha_A*D_A+alphaA*DA-deltaMA*MA && 
MR'== alpha_R*D_R+alphaR*DR-deltaMR*MR && 
A'== betaA*MA+thetaA*D_A+thetaR*D_R 

-A*(gammaA*DA+gammaR*DR+gammaC*R+deltaA) && 
R'== betaR*MR-gammaC*A*R+deltaA*C-deltaR*R && 
C'== gammaC*A*R-deltaA*C 

(b) Uppaal automaton representation. 



Figure 5: Dynamics of genetic oscillator. 



The ODE system can be interpreted as a system behavior near the thermodynamic limit (infinite 
population sizes while maintaining the same concentrations). Alternatively this oscilator can be modeled 
as a system of stochastic chemical reactions, where each molecule is counted as discrete entity moving 
according to Brownian motion laws. Using a standard translation between deterministic and stochastic 
semantics of chemically reacting systems (e.g., Gillespie's algorithm |[T9l ) the coefficients in ODE can 
be interpreted as reaction rates. The reactions are enumerated in Fig. [6al Each reaction is then modeled 
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Figure 6: Stochastic model of genetic oscilator. 



as a separate SHA process shown in Fig. [6bl which can be viewed as an encoding of continuous-time 
Markov chain process. For example, the first reaction means taking one molecule of each A and Da 
substance and producing D'^ with a rate of Ja, which can be interpreted as a transition requiring positive 
amount of A and DA (modeled as integer variables), consuming one molecule each and producing one 
D_A and the reaction rate is gammaA and proportional to available quantities of A and DA. The resulting 
trajectories of quantities A, C and R are displayed in Fig. |7bJ where the patterns seem to resemble the 
one from ODE model (Fig. |7a|, but the inherent stochasticity result in shaky lines (even saw teeth) and 
unpredictably fluctuating amplitudes and their periods. Interestingly, our formalism is flexible enough 
to accommodate a hybrid model combining stochastic aspects as in Fig.|7b]a?i(i continuous aspects as in 
Fig. ITU 

The amplitude of each protein quantity can be measured by the queries E [<=75 ; 2000] (max : q) , 
where 75 is the time limit for simulation, 2000 is the number of simulations and q is either A, C, or R. 
The upper plots of Fig. [8] show the probability density for a range of possible values of amplitude with a 
vertical line for average value. 

Uppaal-SMC can also estimate a distance between peaks by using techniques developed for MITL 
(Metric Interval Temporal Logic - a more expressive property language than a subset of CTL supported 
by Uppaal). The idea of the approach is to translate MITL formula into a monitoring automata which 
start an auxiliary clock x with a first peak and stop with a second peak fTOl. To detect peaks of A when 
its amount rises above 1100 and drops below 1000 within 5 time units, we use the formula (in the tool 
syntax): true U [<=1000] (A>1100 & true U[<=5] A<=1000) . Then the distance between peaks 
can be estimated by measuring maximal value of clock x. The result is shown as logarithm of probability 
density plots in a second row of Fig.[8l The plots show that in most cases the measured distance between 
peaks is about 24.2 hours (slightly more than one day-night cycle). Then there are some smaller bumps 
with several magnitudes lower probability which can be explained by either a) false positive peak as 
MITL monitor is confused by a sudden stochastic saw tooth in signal A, or b) missing a peak or two, or 
even three (in C) if the peak is not high enough to be registered, hence the next one is registered instead. 
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Figure 7: Uppaal-SMC simulations: simulate 1 [<=75] { A, C, R }. 
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Figure 8: Estimated probability density distributions for amplitude and period. 



7 Conclusions 



This paper presents an extensions of the Uppaal-SMC model checker supporting statistical model check- 
ing for stochastic hybrid automata, where dynamic of systems can be specified with ODEs. This is a 
major advance in comparison to existing SMC checkers that can only handle simple derivatives. Cur- 
rently, the tool applies the Euler integration method. In the future, we plan to implement more robust 
methods, such as Runge-Kutta's. Another contribution will be to support rare-events as it is the case in 
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